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FAST  MULTIPOLE  METHOD  IN  SIMULATIONS  OF  AQUEOUS  SYSTEMS 


James  N.  Glosli  and  Michael  R.  Philpott 
IBM  Research  Division,  Almaden  Research  Center 
6S0  Harry  Road,  San  Jose,  CA  95120-6099 


Abstract 

The  fast  multipole  method  (ftnm)  for  calculating  electric  fields  developed  by 
Greengard  and  Rokhlin  (J.  Comp.  Phys.  73,  325(1987)),  has  been  implemented 
specifically  for  molecular  dynamics  simulations  of  electrochemical  problems 
including  boundary  conditions  associated  with  metal  electrodes.  This  order  N 
(number  of  charged  particles)  algorithm,  is  known  to  be  computationally  much 
more  efficient  than  direct  or  Ewald  sum  methods  (order  N2)  for  systems  with 
as  few  as  one  thousand  charged  particles  (equivalent  to  300  water  molecules). 

Timings  and  accuracy  estimates  with  system  size  N  (16  £  N  <  16000)  are  given 
to  illustrate  the  effectiveness  and  efficiency  of  the  fmm. 

As  an  application  of  faun  we  describe  it's  use  with  constant  temperature  mo¬ 
lecular  dynamics  to  calculate  die  dielectric  constant  of  die  spec  water  model  in 
bulk  at  temperatures  298K  and  36  IK.  System  sizes  of  27, 64, 125  and  216  water 
molecules  were  considered.  Comparison  with  the  Ewald  and  reaction  field 
methods  was  made.  At  298K  the  dielectric  constant  was  calculated  to  be 
e=75.5±5%  and  at  361K  £=57.314%.  Both  values  compare  well  with  exper¬ 
iment  and  the  reaction  field  theory  simulations. 

I.  INTRODUCTION 

Molecular  dynamics  and  Monte  Carlo  simulations  of  electrolyte  systems  require  the  evaluation 
of  the  superimposed  Coulomb  fields  of  thousands  particles  (  a  mere  5  nm3  of  water  contains 
approximately  4200  atoms),  fa  this  paper  we  describe  the  faun  of  Greengard  and  Rokhlin  1-4 
and  its  successful  application  to  molecular  dynamics  simulations  of  aqueous  systems  including 
die  calculation  of  electric  polarization  in  model  water  systems.  Specialization  to  interfacial 
problems  of  interest  in  electrochemistry  is  presented  in  outline  only  by  way  of  a  discussion  of 
boundary  conditions  imposed  by  various  interfaces. 

A  number  of  methods  are  in  use  to  calculate  long  range  Coulomb  fields.  Direct  summation  with 
a  cut  off  after  about  1.00  nm  is  common  in  many  commercially  available  codes  like  Charmm  5 
and  Amber  6<  7.  Also  available  are  empirical  recipes  like  the  Hingerty  function  8  and  other 
distance  dependent  dielectric  'constants'  6*  9.  A  brief  survey  of  empirical  dielectric  recipes  in 
use  in  biological  simulations  has  been  given  by  Friedmann  and  Honig  10.  Widely  used  when 
rigor  and  accuracy  and  needed  is  the  Ewald  method  1  1'13.  Also  there  are  a  host  of  plane  wise 
summation  methods  with  origins  in  optics  of  crystals14.  For  homogeneous  systems  the  reaction 
field  method  of  Barker  15  is  unrivaled  for  its  simplicity  and  ease  of  use.  Unfortunately  this 
method  is  not  easily  applied  to  liquid-solid  interfaces  though  attempts  to  extend  it  have  been 
described  161  17.  Ewald  summation  crudely  applied  is  proportional  to  N2  and  at  best  N3/2  where 


N  is  the  number  of  charges.  The  fast  multipole  method  developed  by  Greengard  and  Rokhlin 
1-4  »s  order  N,  and  clearly  the  only  viable  method  for  large  simulations.  The  fmm  technique  is 
also  attractive  for  the  ease  of  implementation  of  a  variety  of  boundary  conditions  such  as  peri¬ 
odic,  Dirichlet,  Neumann  and  mixed  boundaries.  The  method  is  a  recursive  algorithm  based  on 
multipole  expansions  for  the  evaluation  of  the  Coulomb  fields.  The  error  in  the  algorithm  is 
well  controlled  and  can  be  made  arbitrary  small  by  increasing  the  number  of  terms  retained  in 
the  multipole  expansions.  Additionally  an  adaptive  version  of  the  algorithm  is  available  3  in 
which  regions  of  low  or  no  charge  density  are  not  subdivided  when  the  charge  count  falls  below 
a  specified  integer.  This  results  in  considerable  saving  of  computer  time,  and  will  make  the  code 
particularly  useful  in  computations  of  equated  membranes  and  globular  proteins. 

The  method  has  a  number  of  advantages  over  direct  and  Ewald  sum  methods.  As  already 
mentioned  it  is  an  order  N  and  hence  faster  for  large  number  of  particles  (for  a  precision  of 
10*4  the  cross  over  size  varies  from  30  to  1000  depending  on  boundary  conditions).  Also  within 
the  structure  of  the  algorithm  it  is  easy  to  impose  a  variety  of  boundary  conditions.  In  contrast 
to  Ewald  methods  no  addition  computational  complexity  is  introduced  by  imposing  boundary 
conditions.  The  method  is  also  well  suited  to  vector  and  parallel  computation,  allowing  for  the 
possibility  of  exploring  very  large  systems. 

II.  DESCRIPTION  OF  THE  FAST  MULTIPOLE  METHOD 

In  the  calculations  described  later  in  this  paper  the  three  dimensional  version  of  the  fmm  method 
is  used  to  evaluate  the  electrostatic  interaction,  however  for  the  purposes  of  describing  the 
methods  it  suffices  to  consider  the  two  dimensional  version.  The  generalization  to  three  di¬ 
mensions  is  relatively  straight  forward. 

Suppose  simulation  cell  C  contains  N  charged  particles.  The  method  divides  C  (which  we  will 
call  a  level  0  box)  into  four  level  1  boxes.  See  Figure  1.  These  boxes  are  then  subdivided  into 
four  more  boxes  creating  an  array  of  level  2  boxes.  This  process  is  repeated  until  some  pre¬ 
scribed  level  L  is  reached.  In  Figure  1  displays  L=3.  The  finest  level  L  is  chosen  to  minimize 
computation  and  in  doing  so  is  proportional  to  log(N).  The  boxes  formed  by  the  division  of  a 
box  are  called  children  of  the  undivided  'parent'  box.  The  far  field  due  each  to  box  at  the  finest 
level  is  calculated  in  terms  of  it's  multipoles,  up  to  order  n.  The  precision  of  the  method  is  set 
by  the  order  n.  Next  the  multipole  moments  of  the  parents  are  created  from  the  children's  mo¬ 
ments  by  a  linear  transformation,  starting  at  the  finest  level  and  working  upward  to  level  0. 

The  strategy  now  will  be  to  calculate  a  local  expansion  for  each  box  that  represents  the  far  field 
potential.  Beginning  at  level  zero,  the  question  is  asked,  which  boxes  are  well  separated  from 
the  main  computational  cell?  The  answer  will  depend  on  the  boundary  conditions.  For  a  system 
with  a  'free'  boundary  there  are  no  well  separated  level  0  boxes  because  the  there  are  no  other 
level  0  boxes!  For  any  periodic  boundary  condition  it  is  the  array  of  all  images  except  adjacent 
images.  The  local  expansion  for  the  level  zero  box  is  found  by  transforming  the  multipole 
expansion  of  each  of  the  well  separated  boxes.  In  the  case  of  the  free  boundary  this  is  trivial 
since  there  are  no  well  separated  boxes.  However  in  the  case  of  a  periodic  boundary  condition 
there  is  an  infinite  number  of  image  boxes.  Fortunately  in  this  case  symmetry  enables  the  local 
expansion  to  be  calculated  easily.  The  local  expansion  of  a  level  1  box  can  be  considered  to 
consist  of  two  parts.  One  part  representing  particles  included  in  the  local  expansion  of  the  parent 
and  another  part  for  all  particles  in  well  separated  boxes  that  are  not  present  in  the  first  part 
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Figure  1.  Schematic  representation  of  the  L  =  3  subdivision  of 
the  computational  box  C  used  in  the  fast  multipole  method. 

The  number  of  these  well  separated  boxes  is  finite  and  bounded.  The  contribution  to  the  local 
expansion  by  particles  in  the  first  part  is  simply  found  by  transforming  the  origin  of  the  parent's 
local  expansion  to  the  center  of  the  box.  The  second  part  is  calculated  by  transforming  a  finite 
number  of  multipole  expansions.  Once  this  procedure  is  applied  to  all  level  1  boxes  it  is  applied 
recursively  to  level  2  and  so  on.  Finally  a  local  expansion  of  all  particles  in  well  separated  boxes 
is  found  for  each  box  at  the  finest  level.  What  remains  to  be  done  is  to  find  die  potential  due 
to  particles  in  adjacent  boxes  and  this  is  calculated  by  doing  a  direct  sum  over  particles  in  these 
boxes. 

The  Figure  2  schematically  illustrates  the  two  dimensional  analogs  of  systems  of  interest  in 
electrochemistry.  Arrows  show  location  of  metal  boundaries.  In  the  first  (a)  the  simulation  cell 
is  replicated  throughout  2-space,  to  represent  an  infinite  bulk  fluid.  In  cases  (b-d)  it  is  replicated 
laterally  in  the  x  direction,  however  die  y  direction  has  to  be  handled  differendy  in  each  case. 
In  case  b  the  system  is  a  either  an  isolated  Cfree’)  film  or  one  bounded  by  two  low  dielectric 
surfaces.  In  case  c  the  film  is  bounded  by  two  metals  which  require  the  inclusion  of  images 
boxes  C  (generated  from  C  by  charge  conjugation:  q  -» -q  and  reflection  in  the  xy  plane:  z-zQ 
— »  zQ-z  ).  The  bottom  right  case  (d)  represents  an  emersed  electrolyte  film  of  finite  thickness 
on  a  metal  electrode.  In  this  case  there  is  only  one  set  of  image  boxes  C 

The  implementation  of  these  various  boundary  condition  is  straight  forward  in  the  finm  tech¬ 
nique.  For  instance  consider  the  periodic  boundaries.  In  this  case  the  cell  C  is  replicated  in 
space,  as  illustrated  by  Figure  2a.  Since  we  know  the  multipole  expansions  for  the  cell  C,  we 
thereby  know  the  expansions  for  each  of  the  boxes  in  any  image  cell.  To  begin  the  calculation 
of  the  local  expansion  of  cell  C  (i.e.  the  level  0  box)  we  need  to  sum  over  all  transformations 
of  the  coefficients  of  well  separated  level  0  image  boxes.  To  sum  over  all  these  transformations 
each  time  would  be  extremely  costly  in  cpu  time.  However  if  we  use  the  fact  that  the  trans¬ 
formation  operator  is  a  linear  operator  and  the  coefficients  of  each  cell  are  identical  then  we  need 
only  sum  over  all  the  coefficients  to  find  a  new  operator  that  then  can  be  applied  to  the  multipole 
coefficients  of  cell  C  (at  level  0)  in  order  to  find  the  local  expansion  coefficients  of  cell  C.  The 
beauty  of  this  approach  is  that  this  operator  is  independent  of  the  particle  positions  and  therefore 
only  needs  to  be  calculated  once  in  a  lifetime. 
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Figure  2.  Two  dimensional  analogs  of  boundary  conditions  frequently  met  in 
molecular  dynamics  simulations  of  electrochemical  systems.  Arrows  show  lo¬ 
cation  of  metal  boundaries. 


In  die  case  of  an  isolated  film  as  illustrated  in  Figure  2b,  the  same  approach  as  describe  above 
will  work.  The  only  difference  is  that  the  sum  over  transformation  operators  is  restricted  to  a 
single  row.  So  the  summed  operator  will  be  different  from  above  but  still  only  needs  to  be 
calculated  once  in  a  lifetime.  The  introduction  of  metal  boundaries  introduces  images  that  arc 
related  to  the  root  cell  by  charge  conjugation  and  inversion  operations.  For  example  consider 
Figure  2c,  where  a  periodic  film  is  constrained  between  two  metal  plates.  The  images  cells  now 
consists  of  alternating  rows  of  C  and  C  cells.  By  considering  the  heavily  outline  box  as  the 
primitive  cell  in  the  system,  the  system  is  mapped  back  to  a  fully  periodic  system.  One  may 
be  tempted  to  use  the  summed  operator  that  was  calculated  for  the  case  2a,  however  one  needs 
to  be  concerned  with  the  conditional  convergence  of  the  operator  sums.  In  case  2a  the  appro¬ 
priate  order  is  to  sum  over  spherical  shells,  since  one  may  argue  that  the  physical  system  is  to- 
tationally  invariant  even  though  the  replicated  lattice  of  cells  is  cubic  (in  3-space)  or  square  (in 
2-space).  However  in  case  2c  the  isotopic  symmetry  has  been  broken  by  the  metal  plates  and 
the  sums  must  be  performed  in  slab  fashion.  The  case  illustrated  by  Figure  2d  can  be  mapped 
to  case  2b. 

III.  VALIDATION 

In  this  section  the  fmm  method  will  be  compared  with  other  methods  for  evaluating  electrostatic 
interactions.  The  comparisons  are  made  among  eleven  systems  spanning  three  orders  of  mag¬ 
nitude  in  size  N  and  with  three  different  boundary  conditions: 

1)  free  boundaries 

2)  periodic  in  x,y,z 

3)  periodic  in  x  and  y  and  free  the  z  direction. 


For  these  comparisons  as  well  as  the  rest  of  the  paper,  this  implementation  of  fmm  method  re¬ 
tained  the  set  of  spherical  harmonics  (T„,  1 0  <  n  <  7)  in  all  expansions.  This  typically  results  in 
a  relative  accuracy  of  the  order  of  10" 4  The  number  of  levels  was  chosen  to  minimize  the 
computation  time  and  is  dependent  on  the  number  of  particles.  To  a  good  approximation  this 
minimum  is  found  at  L-  log*(A764). 

For  free  boundaries  the  fmm  is  compared  with  direct  summation  methods.  For  the  purpose  of 
this  comparison  it  is  assumed  that  die  direct  method  is  exact  (actually  only  accurate  to  10~13). 
For  the  fully  periodic  system  the  comparison  is  made  with  the  standard  three  dimensional  Ewald 
method.  In  the  two  dimensional  periodic  system  die  two  dimensional  slab  Ewald  method  of 
Rhee  et  al  18  is  used  for  comparison.  In  the  periodic  systems  (two  and  three  dimensional)  the 
accuracy  of  the  Ewald  method  is  controlled  by  the  number  of  terms  in  the  real  and  reciprocal 
space  sums.  Using  a  simple  estimate  of  the  error  the  cutoff  radius  in  real  and  reciprocal  space 
was  chosen  to  achieve  an  accuracy  of  approximately  ID-4.  To  test  the  accuracy  of  both  the  fmm 
and  Ewald  methods  a  high  accuracy  Ewald  method  was  used.  All  the  algorithms  were  coded 
in  double  precision  and  run  on  an  IBM  RS6000  model  540  workstation. 


To  compare  the  methods  an  N  particle  test  system  was  generated  by  choosing  from  a  uniform 
random  distribution  on  the  interval  (-&,V4)  the  charge  and  x,y,z  coordinates  of  each  of  the  par¬ 
ticles.  The  systems  are  made  neutral  by  adding  a  constant  correction  to  each  charge.  In  Table 
1  the  cpu  time  t,  absolute  relative  error  in  the  energy  eE.  and  rms  deviation  of  the  force  eF  are 
listed  for  the  various  methods.  The  cpu  time  t  is  in  seconds  and  eF  is  defined  as. 


eF  = 


(1) 


The  far  right  hand  column  of  table  1  list  the  number  of  levels  L  used  in  the  fmm  calculation. 
At  the  end  of  the  table  is  a  summary  of  a  fit  of  the  execution  time  to  the  form  f =AJV *.  The 
coefficient  p  for  the  fmm  method  is  close  to  one,  as  predicted,  for  all  the  boundary  conditions. 
In  fact  the  deviation  from  one  can  be  traced  to  surface  to  volume  effects.  From  the  fitting  pa¬ 
rameters  a  critical  size  Nr  at  which  the  fast  fmm  method  becomes  most  efficient  is  estimated  and 
it  too  is  listed  at  the  end  of  table  1.  The  estimate  of  Nc=  17  for  the  two  dimensional  slab 
boundary  conditions  is  lower  than  the  actual  cross  over  value  of  30,  but  the  discrepancy  can  be 
traced  to  surface  to  volume  effects  as  well.  As  expected  for  large  systems  the  fmm  is  much 
more  efficient  than  direct  or  Ewald  methods.  In  particular  for  systems  with  finite  geometry  in 
one  or  more  dimensions,  as  in  the  two  dimensional  slab  case  the  fmm  is  clearly  superior  to 
Ewald  methods. 
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eF 
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eF 
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0.0 

0.0 

0.11 
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0.0 

0.0 

0.02 

0.17 

0.0 

0.0 

0.0 

0.0 

0.07 

0.26 

0.0 

0.0 

0.0 

0.0 

0.29 

0.53 

0.0 

0.0 

0.0 

0.0 

1.12 

1.28 

0.0 

lx  10-5 

0.0 

3x  10*5 

4.6 

2.3 

0.0 

1x10-5 

0.0 

3x  10-5 

18.3 

6.2 

0.0 

lx  10-* 

0.0 

3x  10-‘ 

73 

14 

0.0 

3x  IQ*3 

0.0 

lx  10-4 

periodic  3d 


oeriodic  2d 


Ewald 


0.010 


5.9 


lx  10-’ 


8X10-5 


255 


fmm 


0.153 


2xl0-<  4xlO-J 


8x  10-4  2x  10-4 


0.034 


6x  10-4  lx  10-5 


6x  10-4  lx  10-4 


4xlO-J  2xl0-4 


3x  10*4  6x10-5 


2x10-5  4x10-5 


6xl0-4  9x10-5 


7x10-5  3x10*5 


5xl0-4  7x10-5 


2.6 


7x10-5 


2xl0-4  5x10-5 


2xl0-4  4x10-5 


3x  10-4  7x  IQ-5 


10.4 


2x10-5 


6xl0-‘ 


21.5 


lx  10-4  4x  IQ-5 


5x10-5  2x 


Table  1.  The  cpu  time  t  in  seconds  for  an  IBM  workstation  RS6000/540,  relative  error  Ci  in 
energy  and  the  RMS  value  of  the  relative  error  eF  in  the  force. 


N 

■ 

free 

periodic  3d 

periodic  2d 

fl 

Direct 

fmm 

Ewald 

fmm 

Ewald 

fmm 

H 

t 

310 

26 

950 

35 

16000 

<£ 

0.0 

3x  10-5 

2xl0-4 

2xl0-4 

3 

ef 

0.0 

lx  10-4 

2xl0-4 

6xl0-‘ 

Fit 

0 

1.636 
x  10-‘ 

7.234 
x  10-4 

1.839 

xlO-5 

2.905 

xlO-3 

1.077 
x  10-4 

1.774 
x  10-J 

1 

P 

1.960 

1.084 

1.835 

0.976 

1.994 

0.995 

■ 

K 

1069 

393 

19 

1 

IV.  APPLICATION  TO  ELECTROCHEMICAL  SYSTEMS 

Electrochemical  systems  prevade  nature.  At  the  simplest  level  they  can  be  visualized  as  a  col¬ 
lection  of  polar  entities  held  together  by  electrostatic  interactions  and  held  apart  by  the  Pauli 
repulsion  interaction  of  their  constituent  electrons.  Coulomb  interactions  being  essentially  long 
range  in  nature,  control  the  fluctuations  in  electric  polarization  that  manifest  as  the  dielectric 
constant  Since  a  major  component  of  the  total  interaction  energy  is  electrostatic  it  is  essential 
to  compute  this  part  accurately  if  polarization  effects  are  important  In  the  light  of  these  intro¬ 
ductory  comments  it  is  appropriate  that  we  demonstrate  the  power  of  die  fmm  code  as  imple¬ 
mented  by  us,  in  a  calculation  of  the  dielectric  constant  First  however  it  is  necessary  describe 
briefly  the  models  and  interaction  potentials  used  in  these  studies. 

The  sdcc  water  model l9.  In  this  model  water  is  represented  as  a  lennard-jones  atom  centered 
or  the  oxygen  atom  O  with  a  =  0.3160  nm,  e  =  0.6502  kJ/mol.  Embedded  in  the  sphere  are  two 
point  masses  H  carrying  charge  q=0.4238lel  in  electron  units.  The  HOH  bond  angle  is  109.5° 
and  the  OH  bond  length  is  0.1  nm. 

The  total  interaction  energy.  The  Coulomb  interaction  between  molecules  was  represented  as 
sum  of  1/r  interactions  between  all  point  charges.  The  short  range  part  of  the  intermolecular 
interaction  was  modeled  by  a  lennard-jones  potential  between  the  O  atoms  of  each  water,  and 
smoothly  truncated  to  zero  at  R  =  0.72  nm  by  the  function  T  defined  below.  The  electrostatic 
coupling  constant  K  had  the  value  138.936  KJ«nm/(mole*  e1)  in  the  units  used  in  this  calculation. 

"V  complete  interaction  energy  U  is, 


where  i  and  j  are  molecular  indices,  and,  a  and  P  are  atomic  indices.  The  symbol  At  represents 
the  set  of  all  atoms  of  molecule  i  .  The  symbol  R#  is  the  distance  between  the  center  of  mass 
of  molecules  i  and  j.  The  symbol  is  the  distance  between  atoms  a  and  p. 

The  form  of  the  truncation  function  T  is  given  by. 


T(R)  = 


r<rI 

rtl<r<rI 

*Tu<R 


(3) 


where  Rl= 0.68  nm  and  RTU=0.12  nm.  The  intege  m  and  n  control  the  smoothness  of  the  trun¬ 
cation  function  at  R[  and  Rj,  respectively.  In  this  calculation  n  —  m  =  2  which  insured  that  energy 
was  smooth  up  to  and  including  first  spatial  derivatives. 

Bond  lengths  and  angles  were  explicitly  cf  .istrained  by  a  quaternion  formulation  of  the  rigid 
body  equations  of  motion  2^"22.  The  equations  of  motion  were  expressed  as  a  set  of  first  order 
differential  equations  and  a  fourth  order  multi-step  numerical  scheme  was  used  to  integrate  them 
in  time  steps  of  2  fs.  At  each  time  step  a  small  scaling  correction  was  made  to  the  quaternions 
and  velocities  to  correct  global  drift.  Additionaly  the  global  center  of  mass  velocities  in  the  x, 
y  and  z  directions  were  set  to  zero  at  each  time  step  by  shifting  the  molecular  translational  ve¬ 
locities. 


V.  DIELECTRIC  CONSTANT  OF  WATER 

As  a  test  of  the  fast  multipole  method  in  an  actual  molecular  dynamics  simulation  we  used  it 
to  evaluate  the  dielectric  constant  of  water  in  bulk.  There  tw/e  been  numerous  simulations  of 
water  dielectric  constants.  Unfortunately  many  omit  important  details  that  make  it  impossible 
to  repeat  the  calculation.  There  are  however  a  few  that  have  provided  these  details  and  we  have 
selected  one  clear  enough  to  make  comparisons.  Previously  Reddy  and  Berko witz  23  used  the 
reaction  field  method  to  estimate  dielectric  constant  of  spce  water.  In  this  study  we  have  used 
both  fmm  and  Ewald  techniques  to  estimate  dielectric  constant  of  spce  water.  As  necessary  with 
all  simulations  of  Coulomb  fields  the  boundary  conditions  at  infinity  need  to  be  clearly  specified. 
We  have  assumed  that  the  primitive  cell  has  been  periodically  replicated  to  fill  a  sphere  with 
very  large  radius.  Outside  of  the  sphere  space  is  assumed  to  be  filled  with  a  perfect  conductor. 
The  appropriate  fluctuation  expression  for  the  dielectric  constant  of  the  system  with  this 
boundary  conditions  is: 


r_  <M2> 
IVkT 


+  eo 


(4) 


where  M  is  the  total  polarization  of  the  sample  and  e,  is  the  dielectric  constant  of  a  vacuum. 
The  simulations  where  done  for  systems  containing  27,  64,  125  and  216  spce  water  molecules 
at  temperature  of  298K  and  36 IK.  .\11  the  simulations  ran  for  at  least  1  ns  except  for  the  27 
and  64  molecule  systems  at  36 IK  which  were  integrated  0.5  ns.  In  the  simulations  at  the  higher 
temperature  36 IK  we  evaluated  coulomb  fields  by  fmm  and  Ewald  techniques  in  order  to  di¬ 
rectly  compare  the  two  methods.  In  estimating  the  errors  in  the  dielectric  constant  the  temporal 
correlation  between  the  polarization  of  different  samples  is  important  In  this  study  configure- 


tions  were  stored  every  0. 1  ps  and  then  binned.  For  large  bin  sizes  the  samples  should  be  stat¬ 
ically  independent  This  was  tested  by  examining  the  variance  of  the  average  polarziation  as  a 
function  of  bin  size.  If  this  appeared  to  be  approaching  an  asymptotic  value  then  the  bin  size 
was  deemed  to  be  large  enough  ai.J  the  asymptotic  value  of  the  variance  was  used  to  estimate 
the  error  in  t .  Using  this  method  we  found  agreement  between  the  finm  and  Ewald  methods  for 
the  si/.*  of  the  error.  However  our  error  estimate  was  much  larger  that  the  error  estimates  given 
by  Reddy  and  3erkowitz  23  for  runs  of  similar  length  in  their  reaction  field  study. '  However  the 
method  they  used  to  determine  the  error  is  not  clear  from  reading  their  paper. 
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Figure  3.  The  dependence  of  the  relative  dielectric  constant  e/e0  of  spce  Wc.ter 
on  the  number  N  of  molecules  in  the  primative  simulation  cell. 


Figure  3  shows  a  plot  of  dielectric  constant  versus  size  N  the  number  of  spce  water  molecules 
in  the  basic  simulation  cell  for  two  temperatures.  The  different  data  symbols  refer  to  the  results 
from  using  fmm  and  Ewald  methods  to  evaluate  the  electric  forces  acting  on  the  charges.  Oth¬ 
erwise  the  two  codes  were  identical.  The  first  thing  to  note  is  that  the  fmm  and  Ewald  tech¬ 
niques  agree  within  statistical  error.  The  216  molecule  simulation  with  fmm  can  be  compared 
to  Reddy  and  Berkowitz  23  reaction  field  simulation.  Using  fmm  the  value  e=75.5±5%  at 
T=298K  and  e=57.3  ±  4%  at  T=361K  is  found  for  the  dielectric  constant  of  the  spce  water 
model.  In  the  reaction  field  study  e(T=298K)=70.7±l%  and  E(T=373K)=51  ±  .7%.  Using  a 
linear  interpolation  this  would  correspond  to  of  £=54.1  at  36 IK.  At  the  higher  temperature 
T=361K  there  appears  to  be  only  weak  size  dependence,  whereas  at  T=298K  the  size  dependence 
is  clearly  apparent  especially  in  the  range  from  125  to  216  molecules.  The  calculated  values 
agree  well  with  experimental  values  24. 


VII.  CONCLUSIONS  AND  OUTLOOK 

This  paper  described  the  first  electrochemical  application  of  the  fast  multipoie  method  in  a  mo¬ 
lecular  dynamics  calculation  of  the  dielectric  constant.  The  fast  multipoie  method  provides  a 
fast  and  accurate  technique  for  the  evaluation  of  electrostatic  interactions  and  clearly  is  the 
prefered  method  in  systems  with  large  numbers  of  charged  particles.  In  a  system  of  216  spce 
water  molecules,  which  has  648  charges,  the  fmm  took  62%  of  the  time  of  an  Ewald  method 
of  the  same  accuracy.  For  512  spce  waters  fmm  takes  only  30%  of  the  time  for  Ewald.  The 
method  is  well  suited  to  handle  a  variety  of  boundary  conditions.  Given  that  this  method  is 
available  in  adaptive  and  parallelizable  versions  it  is  a  very  powerful  tool  for  use  in  classical 
simulations  of  systems  with  aqueous  phases.  Important  future  applicatons  will  not  be  confined 
just  to  electrolytes  in  the  bulk  or  at  planar  surfaces.  It  is  quite  obvious  that  the  method  can  be 
used  in  simulations  where  water  surrounds  or  is  confined  by  irregular  boundaries  as  are  fre¬ 
quently  encountered  in  biological  systems. 
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